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Abstract 

Imaging Atmospheric Cherenkov Telescopes (lACTs) have resulted in a break- 
through in very-high energy (VHE) gamma-ray astrophysics. While early lACT 
installations faced the problem of detecting any sources at all, current instruments 
are able to see many sources, often over more than two orders of magnitude in 
energy. As instruments and analysis methods have matured, the requirements for 
calibration and modelling of physical and instrumental effects have increased. In 
this article, a set of Monte Carlo simulation tools is described that attempts to 
include all relevant effects for lACTs in great detail but aims to achieve this in 
an efficient and flexible way. These tools were originally developed for the HEGRA 
lACT system and later adapted for the H.E.S.S. experiment. Their inherent flexi- 
bility to describe quite arbitrary lACT systems makes them also an ideal tool for 
evaluating the potential of future installations. It is in use for design studies of CTA 
and other projects. 

Key words: Air showers. Imaging atmospheric Cherenkov technique. Gamma-ray 
astronomy 
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1 Introduction 

Any simulation of the I ACT technique consists of two major components: the 
development of the extensive air showers and emission of Cherenkov light by 
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Fig. 1. Cherenkov light production in CORSIKA for different primary particle ex- 
amples in simulations for a site at 2200 m altitude. Darkness of the particle tracks 
shown increases with increasing emission of Cherenkov light. 



the shower particles as one part and the detection of this light and recording 
of signals by the instrument as the other part. While earlier simulation tools 
used to come with a rather primitive description of the instrumental response 
included with the shower simulation, the simulation tools described here sepa- 
rates these parts into distinct programmes. This allows for enhanced flexibility 
in the implementation and sharing of components in the community. 

The flrst major component is based on the CORSIKA program f]}]. A flrst im- 
plementation of Cherenkov light in CORSIKA, by M. Rozanska, F. Arqueros, 
and S. Martinez, implemented a rectangular grid of detectors - for the non- 
imaging HEGRA AIROBICC instrument [2|. The current implementation, by 
the author of this paper and with contributions by others, unifled the origi- 
nally separate functions for Cherenkov emission by electrons/positrons and by 
other particles, improved its efliciency and accuracy, and added a great level 
of flexibility. The CORSIKA code base includes diflFerent output formats now, 
some more generic and some dedicated to speciflc experiments. Only the most 
generic of these, the lACT option also developed by us, will be covered here. 
With the lACT option, a telescope or other Cherenkov detector in CORSIKA 
is deflned by a flducial sphere containing the reflector or detector. 



The second major component, termed sim_t el array or ^ in its current in- 
carnation for the H.E.S.S. experiment [3] - sim_hessarray, implements all 
the details of the detectors. This includes optical ray-tracing of the photons, 
their registration by photomultiplier tubes, switching of discriminators or com- 
parators at the pixel level and at the telescope trigger level, as well as dig- 
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itization of the resulting signals. Since sim_t el array is highly efficient and, 
in a typical configuration, requires only a small fraction of the CPU time 
needed for the shower simulation, it is usually run several times in parallel, 
reading directly from a common output pipe of the CORSIKA lACT option. 
Each instance might correspond to a different viewing direction, a diflFerent 
atmospheric transmission, or even a completely diflFerent mirror and camera 
definition. 



2 CORSIKA and the lACT option 

2.1 Cherenkov light production in CORSIKA 

Cherenkov light production in CORSIKA (see Figure [1]) with its various op- 
tional] (like lACT, ATMEXT, CEFFIC, CERWLEN, ...) has to be enabled 
before compilation (see for details). This results in a call to the Cherenkov 
subroutine CERENK for each track segment of charged particles as they are 
handled by the CORSIKA particle transport code. The transport code takes 
care of possible interactions, decay, multiple scattering, bending in the geo- 
magnetic field, and ionization energy loss. As far as Cherenkov light production 
is concerned, each straight track segment is defined by start and end point as 
well as by the initial and final energy of the particle (assuming continuous 
energy loss), and its mass and charge. 

Since the track segments can have lengths up to several kilometres at high 
altitudes, the Cherenkov emission subroutine in CORSIKA has to take care 
of changes in the index of refraction of the air as well as of the change in 
velocity of the particle. Track segments in CORSIKA are generally shorter 
when Cherenkov light production is enabled, in particular with the lACT op- 
tion, where typical multiple scattering angles and bending in the geomagnetic 
field should be smaller than the pixel scale in the cameras. Too long track 
segments could result, for example, in too sharp muon rings or even a shift in 

^ CORSIKA options discussed in this paper include: 

ATMEXT: atmospheric extension, providing tabulated atmospheric profiles to all 
CORSIKA variants and refraction to the Cherenkov light. 

CEFFIC: Very simple treatment of atmospheric transmission, mirror reflectivity, 
and quantum efliciency in CORSIKA itself. 

CERENKOV: the master switch for enabling Cherenkov emission. Without lACT 
option, a rectangular detector array is simulated. 

CERWLEN: the index of refraction and, thus, the emission angle depends on the 
wavelength. 

I ACT: the detectors are deflned individually as flducial spheres and the machine- 
independent output data format readable by sim_telarray gets used. 
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Fig. 2. Average lateral distribution of 300-600 nm Cherenkov light in simulations 
of vertical 100 GeV gamma-rays for a site at 2200 m altitude, with different set- 
tings of the EGS step length parameter STEPFC. Too long step lengths result in a 
systematic excess of light inside the light pool (r < 120 m). 

the shower maximum. It has been verified that the reduced step lengths with 
the lACT option (compared to those used in CORSIKA without Cherenkov 
light production and STEPFC=1.0) do not change the amount of Cherenkov 
light or its average lateral distribution in any noticeable way, i.e. well below 1 
percent. 

The STEPFC parameter of CORSIKA could be used to change step lengths 
of gammas and relative to the default settings in EGS [5], which corre- 
spond to STEPFC=1.0. While longer step lengths could be used to improve 
the processing speed, they result in a systematic excess of Cherenkov light by 
about 20% for STEPFC=10 (see Figure [2j). It is reassuring to find that further 
reductions of step lengths (e.g. STEPFC=0.1) will not result in another sys- 
tematic change. The question of proper step lengths is certainly an important 
issue for any shower simulation program - for CORSIKA the built-in default 
appears to be just the right choice. 

While accurate modelling of the showers is perhaps the prime objective in 
CORSIKA, efficiency is important as well. The efficiency aspect is quite im- 
portant in the case of the Cherenkov emission subroutine since most of the 
CPU time is spent there. For many applications, the most important perfor- 
mance gain is achieved by the fact that Cherenkov photons are not simulated 
one by one but in bunches. This concept of bunches was already present in 
the original implementation [2] , with an automatically adapted bunch size as 
a function of the energy of the primary particle, optimized for the AIROBICC 
experiment. 

For imaging telescopes, an automatic selection of bunch sizes is not so easy. 
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Fig. 3. The effect of different bunch sizes on the image obtained with simulated 
lACTs. Ah samples show the same 1.8° x 1.35° part of the field-of-view of a 24 m 
telescope looking at a 62 GeV 7-ray shower from 117 m core distance. The red cross 
in each panel indicates the true direction of the incident 7-ray (for other symbols see 
also FigurefMl). The colour scale at the bottom indicates the intensities found in each 
0.07° pixel, in units of the average intensity from single photo-electrons (p.e.). Top 
row: No CEFFIC option, atmospheric transmission and photon detection handled 
only by sim_t el array; left to right: bunch sizes of 5, 50, and 500 (recommended: 5 to 
10). Bottom row: With CEFFIC option, detection efficiency applied in CORSIKA: 
bunch sizes of 1, 10, and 100 (recommended: <1). Bunch sizes in the left column 
are small enough to show the actual smooth image, in the middle column artificial 
fluctuations are already apparent, and in the right column the artificial fluctuations 
dominated over the true fluctuations in the shower. 

Showers induced by high-energy primaries can be detected from large dis- 
tances while less energetic showers are typically only detected in the region 
of the almost uniformly illuminated central Cherenkov light pool. The opti- 
mum bunch size is highly dependent on the detector configuration and best 
defined by the user in the CORSIKA inputs file. See [4] for instructions. Too 
low values are always safe, but inefficient, while too high values will result in 
artificial image fluctuations, as illustrated in Figured In a typical configura- 
tion with conventional PMTs, a bunch size of 5 to 10 for a wavelength range 
from 250 to 700 nm should be fairly safe without the CEFFIC option. With 
the CEFFIC option, where the detection efficiency is included in CORSIKA 
itself instead of the detector simulation, bunch sizes above one can generally 
not be recommended. 

For a given bunch size, the track segment is sub-divided into steps such that 
in each step one photon bunch is emitted into a random direction on the 
Cherenkov cone. The actual size of each bunch and the Cherenkov cone open- 
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ing angle, defined by cos^ = l/(/?n), depend on the index of refraction n and 
the velocity v = (3c oi the particle at the emission point. The wavelength A of 
the photons is usually unspecified (unless the CEFFIC or CERWLEN options 
are used) and may be randomized by the detector simulation to follow the 
1/A^ distribution in the predefined wavelength range. 

Where even more detailed simulations are desired, a random wavelength can 
be selected in CORSIKA and the wavelength dependence of the index of re- 
fraction (see Figure [4]) taken into account to determine the emission angle 
and the amount of light, i.e. bunch size. This level of detail is not normally 
needed and has a significant impact on performance, and is therefore only 
enabled on request (CERWLEN option). Normally, the index of refraction n 
is assumed independent of wavelength and should correspond to an effective 
wavelength, by default 400 nm. Depending on the configured atmospheric pro- 
file, the refractivity n — 1 is assumed proportional to the density of the air or 
interpolated from a table where the dependence on humidity and temperature 
may be included, typically based on radiosonde data for the site in question. 
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Fig. 4. Wavelength dependence of the index of refraction [g^ and its approximation 
in CORSIKA with the CERWLEN option. The index of refraction as used in COR- 
SIKA typically corresponds to a wavelength of 400 nm. Also shown is the Cherenkov 
emission angle for a very fast particle under normal sea-level conditions. Maximum 
opening angles are 1.37° and 1.34° for 300 and 550 nm, respectively, under the given 
conditions. 
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2.2 The I A CT/A TMO package 



The original Cherenkov light implementation in CORSIKA itself was designed 
to match a rectangular detector array on a horizontal plane. This doesn't allow 
to match most actual or non-horizontal detector layouts. For large zenith- 
angle observations, the concept of horizontally flat detectors will even have 
difficulties with rectangular telescope arrays. 

An extension package for CORSIKA, written in the C language, overcomes 
these limitations. The lACT/ATMO package, together with its corresponding 
interfaces in CORSIKA, serves several purposes: the geometry of an lACT 
array can be matched to actual configurations, use of a machine- and compiler- 
independent data format, use of tabulated atmospheric profiles, and taking 
atmospheric refraction into account. 

An array of telescopes or other detectors is defined by the (x^, yi^Zi) positions 
of the centres and the radii of any number of fiducial spheres (currently 
limited to 1000 in CORSIKA input handling). For a typical lACT the posi- 
tion would best correspond to the intersection of altitude and azimuth axes 
and the radius would be large enough to enclose the whole refiector. No other 
properties of the instrument are required in this part of the simulation. For 
efficient use of the CPU time invested in shower and Cherenkov light sim- 
ulation, the defined array can be re-used multiple times for each shower at 
random displacements with respect to the nominal shower core. By default, a 
circle of given radius i?c,max will be uniformly covered. Depending on the COR- 
SIKA compile-time configuration this can either be a circle in the horizontal 
detection plane as used in CORSIKA without lACT option or in the shower 
plane perpendicular to the direction of the primary particle and projecting 
them back to the detection level. 



Non-uniform coverage is optionally available for implementation by the user. 
This could be used, for example, for generating low energy showers over a 
smaller core distance range than high energy showers. Event weights, which 
are then needed for the reconstruction of spectra, are kept in the provided 
data format. In the more general case, the non-uniform coverage can support 
variance reduction schemes for specific questions, i.e. reducing the variance in 
some resulting number (like the trigger rate, for example) for a given number 
of showers - or reducing the number of showers required for a given variance. 
Usually, preference has to be given to the core positions more likely resulting 
in a trigger, in order to achieve that (thus often called importance sampling 
(3]) - but the optimum distribution depends on the question and detector 
details and cannot be generalised. 

The nominal core position in CORSIKA is along a straight line on the original 
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direction of the primary particle from its starting position at the assumed top 
of the atmosphere - typicaUy 100 to 120 km high. Low-energy charged parti- 
cles, e.g. electrons of a few GeV, suffer substantial bending in the geomagnetic 
field before their first interaction. For electrons detectable with large lACTs 
the deflection can approach i?c,max and would affect detection probabilities 
even under assumption of an isotropic background. The deflection is therefore 
compensated in the lACT/ATMO package. 

Since the possible intersection of each photon path with many spheres would 
be very CPU intense, a rectangular grid at the detection level is set up and each 
sphere is linked to one or more of the grid cells as illustrated in Figure O Only 
spheres linked to the grid cell, where a photon bunch arrives, are checked for 
possible intersections. This scheme is very effective and the CPU time required 
for accumulating photon bunches hitting any detector is generally negligible. 
In order to reduce memory requirements, a thinning scheme can set in when 
too many bunches get collected for a telescope, allowing for high dynamic 
range simulations. In this thinning scheme, similar to thinning of particles in 
CORSIKA, the number of bunches is reduced but the number of photons in 
each of the remaining bunches is increased in compensation. 

The photon bunches intersecting any of the detector spheres in any of the 
randomly displaced array instances are recorded in a machine- and compiler- 
independent, flexible data format, termed eventio^ which was originally de- 
veloped for other purposes 0]. Photon bunches in this format can have 
an unspecified or a specific wavelength, with negative wavelengths indicat- 
ing that detection efficiencies were already applied (CEFFIC option, bunches 
corresponding to photo-electrons). The output can be written to file, with 
optional automatic compression, or it can be piped into another program, ei- 
ther the telescope simulation directly or into an adapter program which in 
turn pipes the output into several different telescope simulations. This way, 
different telescope configurations, different offset angles of the telescopes with 
respect to the source, and different atmospheric transmission models can be 
simulated at the same time. Hundreds of shower simulations, with thousands 
of telescope simulations, can be processed in parallel on a large computer 
cluster without running into any I/O bottlenecks 0] 

Another aspect of the lACT/ATMO package, supporting the CORSIKA op- 
tion ATMEXT, is the possibility to provide tabulated atmospheric profiles. 
These include the density, the atmospheric thickness, and the effective index 
of refraction as a function of altitude. Such a table would typically be derived 

^ With little more than 10 telescope simulations reading CORSIKA data from the 
same network-mounted RAID array of hard disks, the I/O bandwidth would be 
saturated. Piping CORSIKA data directly into sim_telarray, hundreds of com- 
puter nodes may be running 10 telescope simulations each without running into 
I/O limitations. 
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Fig. 5. Definition of grid cells at the detection level to which a detector sphere gets 
linked for detailed inspection of intersections of photon bunches with the sphere. 
The shadow of a sphere is large enough to include all Cherenkov light emitted up to 
10 degrees from the shower direction and intersecting the sphere. For detectors near 
the CORSIKA observation level even most light emitted at larger angles would be 
recorded. Some rays at extremely large angles not relevant for lACTs (but perhaps 
for fluorescence detectors) could go unnoticed for efficiency reasons. 

from radiosonde data for the site in question or from atmospheric models. 
Since the EGS4-based part of CORSIKA cannot interpolate from this table 
but makes use of several atmospheric layers with piece-wise exponential or 
(in the top layer) a linear density profile, the corresponding parameters are 
fitted simultaneously to the atmospheric thickness and density columns at 
program start-up. Where suitable, the table values are interpolated by taking 
advantage of the close to exponential density profile. 

Since systems of lACTs can achieve an accuracy on the level of a few arc sec- 
onds in the positions of point-like sources, atmospheric refraction of Cherenkov 
light has to be taken into account as well [13]. The refraction correction is pro- 
vided by the lACT/ATMO package when the CORSIKA option ATMEXT is 
enabled. It affects the arrival direction, the arrival position, and the arrival 
time at the CORSIKA detection level. These corrections are evaluated once by 
numerical integration and then interpolated as functions of starting altitude 
and zenith angle. Fortunately, the differential refraction, i.e. the dependence 



on wavelength, is not relevant yet for the current systems of lACTs. 

Scattering of Cherenkov light is optionally available, including both Rayleigh 
and Mie scattering. However, this requires preparation of special data files, is 
very CPU-intense, and thus not normally used. In addition, scattered Cherenkov 
light is irrelevant within the short integration times of lACTs (lo| . The scat- 
tering option is therefore not included in the package version distributed with 
CORSIKA but available on request only. 



3 The sim_telarray telescope simulation package 



3. 1 Overview 



The sim_t el array package was originally developed for the HEGRA I ACT 
system j9|. In this early implementation, many aspects were still hard- wired, 
e.g. that the telescope trigger would be a next-neighbour trigger with a pixel 
logic corresponding to discriminators with immediate switching and sharp 
rising and falling edges of the signal. For the design issues with H.E.S.S. and 
future telescope systems many aspects were eventually generalized and are 
fully configurable at run-time. Each telescope can be separately configured, 
including its optical layout and properties as well as the camera configuration, 
how telescope and system trigger conditions are formed, and also how the 
signals get digitized and recorded. In its current implementation, it is also 
referred to as the sim_hess array package although it is in no way specific 
to H.E.S.S. By exchanging a number of configuration files, any other I ACT 
system can be simulated just as well. 

Realistic and detailed detector simulations are key features of the sim_telarray 
package but efficiency is just as important. Often, simulation studies require 
a billion and more events, sometimes with many telescopes in each event - in 



some studies for CTA pjl close to 100. Full simulation of night-sky background 
(NSB) and electronic noise in every pixel of every telescope event would be 
prohibitively slow for such studies. In other cases, e.g. for finding reasonable 
trigger conditions where random triggers due to night-sky background are at 
a manageable level, the full simulation of each and every event is unavoidable. 

As long as night-sky induced triggers are negligible and most simulated events 
will not trigger a telescope, sim_telarray offers several short-cuts in the sim- 
ulation which can boost efficiency by large amounts without losing events 
that would eventually trigger a telescope. These short-cuts require initial full 
simulations to determine safe lower limits for the numbers of Cherenkov pho- 
tons and the resulting number of photo-electrons needed to trigger a telescope 
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(except for NSB-induced triggers). Where the required number of photons ^ 
based on the sum of aU photon bunches hitting the fiducial sphere - is not met, 
the simulation can be bypassed completely. Otherwise the optical ray-tracing 
and photon detection are carried out and the sum of all photo-electrons is 
compared to the second limit. Only where that is met as well we have to go 
through the detailed simulation of the electronics response. The difference be- 
tween complete bypassing and full simulation without much Cherenkov light 
can be up to 5 orders of magnitude in time, with little more than 0.01 seconds 
for the full simulation of a H.E.S.S.-I type telescope on a fast CPU core. In 
real productions, the short-cuts typically result in improvements by factors of 
10 to 100. 



3.2 Optics simulation 



All lACTs to date have segmented reflectors made with mirror tiles of spherical 
type. The optical ray-tracing in sim_telarray fully reflects this current usage, 
with spherical mirror tiles of round, hexagonal, or square shape. The optical 
quality of the tiles can be adjusted to laboratory measurement, including 
variations in focal length or non-perfect surfaces. Also the mirror adjustment, 
typically performed with star light falling onto the lid of the camera, will have 
a limited accuracy and is taken into account when every individual mirror 
tile is conflgured at program start-up. How well this approach can represent 
real telescopes is demonstrated in the case of the H.E.S.S. I reflectors, based 
entirely on mirror quality lab measurements and the measured accuracy in the 



mirror alignment procedure p^p^, without any free parameters. Not only are 
the general characteristics of the optical PSF reproduced but also the detailed 
shapes of star light imaged onto the camera lid [14] (see Figure [6]), at least for 
the zenith angle range where the alignment was carried out. The zenith- angle 
dependence of the alignment quality, due to dish deformations, can be set as 
well in order to match measurements or separate flnite elements calculations. 

The optics is deflned by the overall focal length /tei, a table of mirror positions 
(x^. Hi) perpendicular to the optical axis, shapes, and sizes as well optional Zi 
positions and mirror tile focal lengths If the Zi positions are omitted, mirrors 
will be positioned following variants of either a Davies-Cotton or parabolic 
dish shape, taking mechanical accuracies into account. Mirror tiles can have 
either all the same focal lengths, speciflc fi from the table, or automatically 
optimized fi derived from the desired dish shape and {xi^yi). 

In the case of a parabolic dish, the optimal mirror tiles would clearly be 
paraboloid segments, while all I ACTS so far have been equipped with spheri- 
cal tiles. The optimum focal lengths of spherical tiles on a parabolic dish can 
be obtained by calculating the principal (maximum and minimum) radii of 
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Fig. 6. Measured (left) and simulated (right) distribution of star light imaged onto 
the lid of a H.E.S.S. camera at an off-axis angle of 2.3°, adapted from The 
fields shown are about 10.5 cm (0.40°) wide. 

curvature of the paraboloid at the position of the mirror segment. Half the 
geometric mean, i.e. 0.5 (rmm^max)'^^^, is used for automatically optimized fi 
but the simple distance of the mirror tile to the system focus would be rea- 
sonable as well, as the two lead to almost the same result. Practical aspects 
of the mirror production for a real lACT may dictate to build mirror tiles of 
only a few different focal lengths or even a single focal length, within specifi- 
cations. The optical ray-tracing of sim_t el array is well suited for comparing 
the optical point-spread functions (PSF) of different designs and optimize the 
design. 

An example study for the H.E.S.S. -2 telescope now under construction is il- 
lustrated in Figure [71 demonstrating that the optical PSF is only marginally 
deteriorated by choosing a single focal length for all mirror tiles. 

Since the propagation time of light is included in the optics simulation, the 
time spread due to non-isochronous mirror configurations is automatically 
reproduced. While a parabolic dish guarantees the shortest possible Cherenkov 
pulse durations, different dish shapes may result in a still acceptable time 
spread while improving the off-axis PSF. The sim_telarray program therefore 
allows for variations of the assumed dish shape with respect to either a Davies- 
Cotton or parabolic base design. Since the most important aspect here is 
the radius of curvature of the dish this is largely equivalent also to elliptical 
intermediate dish shapes. 

Mirror tiles can be of hexagonal, square, or circular shape. Since the test for 
intersection with every mirror tile could be very time consuming in the simula- 
tions, a grid scheme somewhat similar to that in the CORSIKA lACT option 
is applied, this time along the dish onto which the mirror tiles are mounted. A 
mirror may be associated with one or more grid cells. By intersecting the path 
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Fig. 7. Study of the point spread function of (Cherenkov) light from a point source at 
10 km distance and 1.5° off-axis imaged to the pixel entrance plane in a H.E.S.S.-2 
type telescope (/ = 36 m, coordinates in units of centimetres). The worst case 
off-axis orientation, along the larger vertical dish diameter, is assumed. Top: ideal 
mirror alignment and individually adapted mirror tiles. Middle: all mirror tiles of 
identical focal length. Bottom: in addition with realistic mirror alignment accuracy. 
R.m.s. widths of the light distributions along radial and transverse directions are 
indicated in the figures, as well as 50%, 80%, and 95% containment radii. Note the 
extremely fine pixels of 0.067°. 



of the incoming photon with the dish one grid cell gets selected and the de- 
tailed intersection has to be evaluated only for the few mirror tiles associated 
with that cell. Every mirror tile is once set up with random misalignments 
and misplacements according to realistic specifications. Its actual focal length 
may randomly vary around the defined or optimized value. The surface is as- 
sumed to exhibit some random roughness resulting in a small smearing of the 
reflected rays. 
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As an optional feature, not normally used for performance reasons, the support 
structure for the camera can be included in the ray-tracing, approximated by 
an arbitrary number of cylinders. This way, the shadowing by the camera 
support structure and by the camera itself and its open lid can be evaluated 
in detail, taking both incoming and reflected photons into account [l2|. For 
normal production, the shadowing by the support structure is parametrized 
as a function of the angle w.r.t. the optical axis and only the shadowing by 
the camera is explicitly included. 

In order to make such optical studies possible, and also for realistic tests of 
the impact of a non-uniform NSB on the shower reconstruction, sim_telarray 
allows to set up additional light sources, normally at inflnite distance (stars) 
but optionally at flnite distance to evaluate the imaging of speciflc positions 
along a shower axis. For the optical studies the listing of reflected photons can 
be written to a flle. For the shower simulations, the additional stars result in 
an enhanced noise in the illuminated pixels, for efliciency reasons determined 
once at program start-up. The average or direct current (DC) of stars in the 
photo tubes - which gets decoupled by high-pass fllters in real detectors ^ is 
subtracted and not present in the signals. 

With mirror tiles aligned to focus star light on the closed lid of the camera, 
the Cherenkov light from the showers is focused behind the lid, at an offset 



d = if-' - D-'r' - f 

for light emitted at a distance D from the telescope of focal length /. See (IB] 
on possible strategies for optimal focusing on showers. 

As a flnal step of the optics simulation, the angular acceptance of the pixels 
is taken into account since pixels are typically equipped by Winston cones or 
such to avoid gaps between neighbouring pixels as far as possible. For practical 
reasons, this is made part of the photon detection efliciency. As far as the ray- 
tracing is concerned, it is suflicient to flnd out which (if any) pixel is hit 
at the entrance of its Winston cone or similar and optionally if the photo- 
cathode at the base of the Winston cone is hit directly or not. The shapes 
of pixel entrance and unobstructed photo-cathode can be hexagonal, square, 
or circular. A more detailed set-up, activated in the camera conflguration 
flle, may include angular acceptance tables, averaged over the pixel area, as 
obtained from measurements or detailed ray-tracing of the light collectors. 
Whenever the relevant data is available, the more detailed mode is usually 
preferred. 
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3.3 Photon detection 

The photon detection probabihty e basicaUy depends on the atmospheric 
transmission T, shadowing factor S by the camera and its support structure, 
the mirror reflectivity i?, the efliciency E of the hght concentrator (optionaUy 
whether a photon hits a pixel directly onto the photo-cathode or via reflection 
on the light concentrator), on the quantum efliciency Q of the photo-cathode 
and flnally the collection efliciency C inside the PMT. Most of them are func- 
tions of wavelength: 

e{\)=T{\)SR{\)E{\)Q{\)C{\). 

While fluctuations in the reflectivity of individual mirror tiles are smoothed 
out, the fluctuations in quantum efliciency have to be taken into account. Such 
fluctuations in quantum efliciency from PMT to PMT are usually obtained 
by measuring the photocathode currents of all PMTs induced by a calibrated 
light source, e.g. the CB (Corning Blue) value provided by the manufacturer. 
The collection efliciency of typically 80 to 90% is the probability that a photo- 
electron actually hits the flrst dynode of a PMT and is effectively multiplied 
rather than elastically scattered. It can be accounted for in two ways. Either 
by including it in the detection probability as above and the pulse amplitude 
distribution at the anode only accounting for photo-electrons multiplied at 
the flrst dynode. Or by including the non-amplifled photo-electrons in the 
amplitude distribution of single photo-electrons {single p.e.)^ see Figure[8l For 
that purpose, the single-p.e. amplitude is randomly drawn according to a table 
rather than from a normal distribution. The small wavelength dependence of 
the collection efliciency is a consequence of the photo-electric effect. It may 
be accounted for by adjusting the quantum efficiency according to 

Q{X)=Q{X)C{X)/C, 

where C is the average or effective collection efficiency. In Figure [9] this is 
illustrated for PMTs evaluated for H.E.S.S. 

For a realistic detector response, both the quantum efficiency and the volt- 
age required to achieve the desired gain can vary from PMT to PMT. For 
simplicity, the same scaling factor is applied to the quantum efficiency at all 
wavelengths. Variations of the voltage U result in variations of the PMT tran- 
sit time which is proportional to 1/VU. The simulation also includes a transit 
time jitter for individual photo-electrons. 

The pixels are permanently exposed to NSB at a level configurable for every 
pixel (including light of additional bright stars) , in units of photo-electrons per 
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Single p.e. amplitude [ peak ampl. ] 



Fig. 8. Single-p.e. amplitude distribution at the PMT anode for the two ways to 
take a 85% PMT collection efficiency into account. Solid line: collection efficiency 
included in the amplitude distribution (the mean amplitude is 85% of the ampli- 
tude in the peak of the distribution). The distribution here is based on a detailed 
simulation of the multiplication process in PMTs with parameters fitted to match 
low-noise measurements between 0.1 and 3 times the peak amplitude. Dashed line: 
collection efficiency is included with the quantum efficiency and a normal distribu- 
tion is assumed. 

nanosecond. For existing detectors, this number is obtained from the image 
noise. An utility program is provided for evaluating the impact of different 
quantum efficiency curves, different mirror reflectivity etc. on the NSB rate, 
assuming the NSB spectrum measured by Benn and Ellison [13] . In real detec- 
tors, the DC (direct current) contributions get decoupled by high-pass filters. 
In the simulation, the mean amplitude of NSB signals is therefore subtracted 
from the signal baseline. For every random NSB photo-electron, the corre- 
sponding pulse shape with randomized amplitude is added on top of that 
baseline. 

Another important NSB aspect taken into account is the afterpulsing typically 
induced by ions in the PMT hitting the photo-cathode. These afterpulses 
typically come hundreds of nanoseconds after the original electron cascade 
by which the atom was ionized. For Cherenkov photons, the afterpulses will 
therefore always be well outside the readout window. NSB photons registered 
hundreds of nanoseconds earlier though can be well responsible for afterpulses 
within the time window relevant for the trigger and signal readout. The ampli- 
tude distribution of afterpulses extends to much larger amplitudes than that 
of normal signals (see Figure [IHD . This fact is taken into account by having 
different single-p.e. amplitude distributions for Cherenkov light and for NSB, 
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Fig. 9. Quantum efficiency (upper curve), quantum efficiency times collection effi- 
ciency (lower curve), and the latter divided by the mean collection efficiency of 85% 
(middle curve) as evaluated for a sample of PMTs for the H.E.S.S. experiment [l^. 
The sim_telarray program can be either configured to include the mean collection 
efficiency in the single-p.e. amplitude distribution (together with the middle curve 
for the effective quantum efficiency) or to include it with the quantum efficiency 
(using the lower curve). 

with the afterpulse amplitude distribution added for the NSB signals only. 
How important proper accounting for afterpulses can be is demonstrated in 
Fig.ini 




2 4 6 8 10 

Single p.e. amplitude [ mean prompt ampl. ] 

Fig. 10. Amplitude distribution for prompt single-p.e. signals and after including 
the afterpulses (fits to two types of PMTs, as compared in Figure [TT]). 
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Pixel threshold [ mV ] 

Fig. 11. Example evaluation of PMT specifications for use in a future Cherenkov 
telescope. One PMT type has a slightly higher quantum efficiency but also a five 
times higher afterpulse ratio than the other (see Figure [TOl). Simulations assume a 
mean single-p.e. amplitude of 23 mV for both and the same night-sky background 
resulting in random camera triggers (without any Cherenkov light). For the thin 
dotted and dash-dotted lines the afterpulsing is ignored, resulting in only a small 
increase in pixel comparator thresholds required with high Q.E. PMTs to achieve 
the same trigger rate as with low Q.E. PMTs. When afterpulsing is included (thick 
solid and dashed lines), the higher Q.E. PMT turns out to be unusable, at least 
with 3-fold pixel coincidences, because the required threshold would be at least 
three times as large as for the low Q.E. PMT type. 

3.4 Electronic signals and trigger 

After the photons are ray-traced to the PMT pixels and there may be a chance 
of triggering of the telescope, based on the number of photo-electrons, detailed 
simulation of the electronics response is started. Several different pulse shapes 
are involved. These begin with the pulse shape for a single p.e. at the input 
of a comparator or discriminator for each pixel. Similar but different pulse 
shapes apply for each channel in the digitization process. Two ways of having 
more than one channel per pixel are supported here: different pre-amplifier 
gains for an enhanced dynamic range (as used with H.E.S.S.) and multiple, 
phase-shifted FADC modules for achieving a high sampling rate (as used with 
HEGRA). 

For the analogue signal at the comparator or discriminator input, an internal 
sampling short compared to actual pulse durations is required. Typically, a 
250 ps time interval is used for pulses of a few nanoseconds. Where the camera 
digital electronics should provide only an integrated signal, a sufficiently fine- 
grained internal sampling is used in the first place and the signal is added up 
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afterwards. Otherwise the actual samphng of the read-out system is used and 
the common phase of the digitization (across the whole camera) is randomized 
with respect to the photon arrival times. For each photo-electron, the pulse 
shapes are scaled by the random single-p.e. response and shifted according to 
the photon arrival time plus PMT transit time including a random jitter. All 
single-p.e. signals from the Cherenkov light as well as the NSB photons are 
added up. For pixels exposed to too much star light, limits are available when 
the pixels are disabled in the trigger and when their HV would be turned off 
to avoid damage. 

For a realistic response of the trigger system, the switching behaviour of the 
comparators or discriminators and the rise and fall times of their output signal 
- typically on the order of a nanosecond - has to be taken into account, as 
illustrated in Figure [121 And they will not switch instantly when the signal 
exceeds the predefined threshold. This can be either achieved by demanding a 
minimum time-over-threshold or a switching charge (minimum area between 
pulse shape and threshold level). The sim_telarray program even takes pos- 
sible variations in the output amplitude from device to device - 10% being 
not uncommon - into account. Since the telescope trigger decision is based (at 
least effectively) on a second comparator or discriminator decision applied to 
the sum of pixel logic outputs, this has testable consequences. A fully digital 
trigger decision would correspond to zero rise and fall times as well as iden- 
tical amplitudes of all pixel logic outputs. In reality, the signals are so short 
that the pulse shape of pixel logic outputs should not be ignored. As an ex- 
ample, the - originally unexpected - smooth change in trigger rate of the first 
H.E.S.S. telescope as a function of the multiplicity threshold could be readily 
explained by the available simulations, with comparator properties according 
to manufacturer specifications. This is illustrated in Figure [T3l 




Fig. 12. Three comparator or discriminator output signals partially overlapping in 
time may or may not be sufficient for a telescope trigger set to a 3-fold multiplicity 
threshold. Detailed treatment of trigger switching is required. 

The telescope trigger logic, i.e. comparator or discriminator outputs from 
which pixels can be combined to form a telescope trigger, is different for each 
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Fig. 13. Telescope trigger rate as a function of the threshold level corresponding to 
the pixel multiplicity. Data from the first H.E.S.S. telescope [18], taken in different 
nights and corrected for zenith angle dependence. Simulations including only proton 
showers were scaled up in absolute rate by a common factor to match the data (no 
absolute calibration being available at that time). The telescope triggers when the 
(analogue) sum of the pixel comparator outputs in any sector of 64 pixels can switch 
the sector comparator with programmable threshold. That threshold is shown here 
in units of the average pixel comparator output, i.e. effectively the number of pixels 
fired at the same time. Idealized simulations would show steps at integer multiplicity 
values, while the realistic simulations show only a very modest flattening between 
integer values. Three sets of simulations were carried out, with random variations of 
quantum efficiencies, comparator output levels, etc., as indicated by manufacturer 
specifications. 

lACT array. In HEGRA a next-neighbour logic was used and in H.E.S.S. a 
sector logic based on 8 x 8 pixel sectors, while Smart Pixel [l^ and other 
designs differ again. The sim.telarray program aims to be fully flexible in 
this area and its configuration files include lists of all pixel combinations from 
which a telescope trigger can be formed via a (multiplicity) threshold in the 
sum of corresponding comparator/discriminator outputs. 

The last step in the trigger decision of an lACT array will be the system 
trigger, usually as a requirement of at least two telescope triggers within a 
short time window, after correcting for the expected time delay as a function 
of the system viewing direction and the telescope positions. The required width 
of the window - typically several 10 ns - depends not only on the fields-of-view 
of the telescopes but also on the stability of the transmission lines (cable or 
optical fibres) from the telescopes to the central or distributed system trigger 
logic. Since random system trigger rates are typically very low and the width 
of the system trigger time window therefore usually chosen wide enough, any 
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changes in transmission delays can be ignored by sim.telarray. Since the 
delay times are compensated automatically in sim_telarray, the telescope 
multiplicity and the width of the coincidence time window are usually the 
only parameters needed. Optionally, sim.telarray allows read-out of non- 
triggered telescopes, as it happened to be the case for the HEGRA lACT 
array. 

3.5 Telescope raw data and other output data 

The main output data from sim_t el array should correspond to what would 
be recorded from real telescopes {raw data) as closely as possible. This can be 
full sample-mode digitized pulse shapes for every pixel or the signal sum over 
the given integration window. Different types of zero-suppression are possible, 
like suppression of low-gain channels when the amplitude in corresponding 
high-gain channels is small or full suppression of pixels with small amplitudes. 

The underlying data format is based on the eventio machine-independent for- 
mat, like the CORSIKA I ACT output itself but using more and different data 
block types. Strictly simulation-related data, like the nature, direction, and 
true energy of the primary particle or the actual displacements of telescope 
arrays with respect to the nominal core position, are kept separate from the 
raw data, i.e. in separate data blocks but in the same data stream or file, and 
can be optionally omitted or filtered out later for realistic data challenges. The 
data format is largely independent of that used by past and current lACT ar- 
rays like HEGRA and H.E.S.S. but conversion to experiment-specific formats 
is generally straight-forward and provided through utility programs. 

Data that is needed for calibration, like conversion from ADC counts to photo- 
electron numbers, laser/LED flat-fielding, muon-ring efficiency etc. can be 
derived from full simulations of the various types of calibration runs. The 
more basic numbers like pedestals, photo-electron conversion factors, and fiat- 
fielding coefficients are also directly provided in the output data. For the 
pedestals, which may depend on a varying temperature and other things in a 
real detector, a random 'measurement' error with configurable r.m.s. is applied 
to the reported value. 

3. 6 Built-in shower reconstruction and event visualization 

The sim_t el array program includes a basic geometrical shower reconstruc- 
tion, based on second moments (Hillas parameters) and stereoscopic recon- 
struction methods originally developed for HEGRA. The shower direction is 
defined by the weighted average of the pair-wise intersection points of Hillas 
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second- moments major axes in every two-telescope combination with two well- 
defined images. Only those reconstruction steps are carried out which do not 
require any look-up of other simulation data. Energy lookup and gamma- 
hadron-discrimination are thus precluded here. Instead, the reconstruction is 
oriented towards immediate visualisation of the data. A dedicated Postscript- 
generator for very compact plot output was developed for that purpose. See 
Figure [HI for an example. 




Fig. 14. Example of the sim_telarray event visualisation, showing the same shower 
(a 7-shower of 460 GeV energy at 190 m core distance) in telescopes with different 
pixel sizes of 0.07 to 0.28°, in a 9-telescope array evaluated for the CTA project [llj]. 
For clarity, only part of the field-of-view of one and the same of the nine telescopes 
is shown for each pixel size. The second moments ellipses are shown at one and 
two times their actual size, with the major axes drawn up to four times their actual 
length. Reconstructed and actual shower directions are indicated by the small circles 
and crosses near the camera centres, respectively. 

4 Conclusions and outlook 

The CORSIKA program with the lACT/ATMO extensions aims to provide 
the most detailed simulation of air showers and the production of atmospheric 
Cherenkov light, both for imaging and non-imaging detectors. For efficiency 
reasons not every possible detail is normally used, like the wavelength de- 
pendence of the index of refraction, but is available on demand. A variety of 
options are at hand to adapt the level of detail and performance as needed 
for specific applications. Parts of the photon detection probability can be ap- 
plied within CORSIKA or can be provided by the second simulation stage, 
the detector simulation. The machine- and compiler-independent eventio data 
format allows for a variety of schemes. The lACT extension, normally but not 
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necessarily using this flexible data format, allows treatment of the largest 
telescope or detector arrays envisaged for coming generations of atmospheric 
Cherenkov experiments or observatories. It is part of the standard CORSIKA 
distribution for several years now. 



For the second stage of the simulation process, a dedicated program for imag- 
ing telescopes with tessellated primary mirrors, termed sim_telarray, aims 
to provide the greatest level of detail in the simulation of all optical and elec- 
tronics components involved in the measurement process. At the same time, 
sim.telarray is highly eflicient. The eventio data from the CORSIKA lACT 
option can be piped concurrently into several sim_telarray processes, each 
with a different conflguration, while avoiding I/O bottlenecks on large com- 
puter clusters. 



The simulation results from sim_telarray can be either written directly into 
an experiment-speciflc format, as initially done with the HEGRA I ACT sys- 
tem, or into a more generic, eventio-hdised format. The resulting data can be 
analysed directly from the latter format or converted into experiment-speciflc 
formats, e.g. the H.E.S.S. data format, and then analysed with the same soft- 
ware as measured data. 



A wide range of cross checks with measured data - mainly by H.E.S.S. - con- 
flrms that CORSIKA and sim_telarray can reproduce the data extremely 
well, for newly commissioned telescopes just based on pre-production lab mea- 
surements of mirror qualities, quantum efliciencies and so on. Examples of that 
include the optical PSF [13], the intensity, ring radius, and ring width of muon 
ring images [20|, image shapes of gamma and background showers as well as 



the reconstructed shower PSF for point-like gamma-ray sources j2lLl22[|. as well 



as the trigger rates [18|]. After commissioning, H.E.S.S. simulations only had 



to be adapted to the slowly degrading optical throughput of the telescopes. 



Future developments in sim_t el array may include non-spherical mirror tiles 
or secondary mirrors, and extensions to the photon detection code, as detectors 
with novel properties come into operation. Since the aim is to provide data like 
a real instrument, more capabilities of future read-out systems, like advanced 
post-processing in FPGAs which may include peak flnding and shape analysis, 
are among the most recent additions. 



The CORSIKA program together with the lACT/ATMO package is available 
from Forschungszentrum Karlsruhe. For details see 
http : //www-ik . f zk . de/corsika/. 



The sim_t el array program is available on request from the author. A release 
under an open source licence is foreseen. 
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